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Abstract 

Pseudospectral collocation methods and finite difference methods have been used for 
approximating an important family of soliton like solutions of the mKdV equation. 
These solutions present a structural instability which make difficult to approximate 
their evolution in long time intervals with enough accuracy. The standard numerical 
methods do not guarantee the convergence to the proper solution of the initial value 
problem and often fail by approaching solutions associated to different initial condi- 
tions. In this frame the numerical schemes that preserve the discrete invariants related 
to some conservation laws of this equation produce better results than the methods 
which only take care of a high consistency order. Pseudospectral spatial discretization 
appear as the most robust of the numerical methods, but finite difference schemes are 
useful in order to analyze the rule played by the conservation of the invariants in the 
convergence. 
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1 Introduction 

In this paper we study from a numerical point of view the so called breather 
solutions of the focusing modified Korteweg-de Vries equation: 

d t u + d*u + 2d x (v?) = 0, i£l, m 
u(x, 0) = uo(x). 

This equation is a canonical non- linear dispersive equation pQ, and therefore it 
appears as a good approximation of different physical problems as the motion 
of the curvature of some geometric flux Ref. ([2], [3], [4]), the vortex patch, 
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ferromagnetic vortices Ref. [SJ, fluid mechanics Ref. ([B], [7J, [5], [5]), traffic 
models, anharmonic lattices, hyperbolic surfaces, etc. 

Travelling wave solutions that exhibit one isolated hump are easily found 
either by direct integration of the corresponding o.d.e. or using the inverse 
scattering method. In the case of the real line they are explicitly given by 

u(x, t) = Bsech{B{x - B 2 t)), (2) 

with B € R. Notice that the above expression could admit another two parame- 
ters to consider the traslations in space and time. As we can see, the travelling 
wave propagates to the right with speed /3 2 . From the inverse scattering point 
of view these solutions are characterized by the property that the reflexion co- 
efficient has a single pole located at the imaginary axis. Solutions of more than 
one hump can also be constructed and they correspond to a reflexion coefficient 
with more than one pole in the imaginary axis. When these poles are differ- 
ent, they represent humps (up or down) of different heights. Therefore, these 
humps travel at different speed and they collide in an almost elastic way, see for 
example [10]. It is particularly interesting the degenerate case when we have 
only one pole, which is double. In this paper we will pay considerable attention 
to this situation considering solutions that evolve asymptotically in time as two 
equal humps such that one is up and the other is down (Figure [Tb). It's widely 
used the term double pole solution to describe such solutions and we shall use 
it in these pages. The numerical simulations that start with this family of ini- 
tial conditions converge to solutions of different type depending on the method 
chosen, and not all of them are satisfactory. It is of fundamental interest to 
clarify why some conventional numerical methods turn an initial condition from 
a family of solutions into a different one and which improvements in the method 
might prevent these instabilities. 

Another relevant family of solutions is the one formed by the so called 
breathers (see Figure QJi). They were firstly obtained by M. Wadati in [TT] 
and describe oscillating pulses that do not disperse. They are determined, up 
to translation in time and space, by two real parameters (a,/3), which corre- 
spond to the frequency of the pulse and the amplitude-width of the envelope. 
The phase velocity of the pulse is given by 3/3 2 — a 2 , and the group velocity by 
fP — 3a 2 , a > 0, so that for a large with respect to j3 the pulse propagates to 
the left with a velocity 3a 2 . In fact, in this case they can be approximated by 

B 2 

u(x, t) w -2— sin[a(a; - (3B 2 - a 2 )t)}sech\B(x - (B 2 - 3a 2 )i)l. (3) 
a 

Wadati's approach to construct these solutions is based on the inverse scat- 
tering method. The breathers are characterized by the fact that the poles of 
the reflexion coefficient, denoted as ±a + i/3, arc symmetric with respect to the 
imaginary axes. This family of solutions for a large a with respect to B was used 
by Kenig, Ponce and Vega in [12] to prove discontinuity of the flowmaps associ- 
ated to mKdV equation in the Sobolev spaces H s of functions with s derivatives 
in L 2 (R). This lack of continuity comes from the fact that two breathers with 
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different speeds can be very close at time zero in a Sobolev norm but, because 
they don't disperse, they will eventually separate and therefore the difference of 
their Sobolev norms becomes very big. The point is that the time of separation 
can be made arbitrary small by taking a large enough. The threshold for this 
lack of continuity occurs for s — 1/4 that turns out to be sharp. 




Fig. 1: (a) Breather type solution of mKdV at time t=0 for a = 7 and /3 = 1. 
(b) Double pole solution for (3 — 1 at time: t — (solid line), t = 50 
(dashed line) and t = 5000 (doted line), after translation to the axis 
origin to show the logarithmic spread in time. 

As far as we know, these solutions exist for the mKdV equation but not for 
the classical KdV equation (i.e. the nonlinearity u x is changed into uu x ). This 
fact gives us a reason to study mKdV better than KdV. 

A natural question that comes up from the observations we have just made, 
is what are the stability properties of these breather solutions when a is large. 
Due to the high oscillations of the pulse, numerical methods based on finite 
difference schemes do not look appropriate, and we will see later on that this 
is the case. However the Fourier transform of the pulse is highly concentrated 
around the frequency a, and therefore pseudospectral methods seem much more 
natural in this case (Ref. [13], [14] and [II])- We will see in the next sections 
that they are in fact very robust. 

One could wonder what happens when finite difference methods are used for 
small a. Notice that the solution is real analytic, and that the available well- 
posedness existence theory for the Cauchy problem (see for example [16 ) allows 
us to conclude that the regularity is preserved along the flow and in principle, 
one shouldn't expect any numerical instability. Nevertheless the construction 
of the breather solution made by the inverse scattering method suggests that 
some instability can arise when we get close to the degenerate case a = 0. The 
reason is obvious. On one hand the double pole solution can be obtained from 
the two soliton solution when the poles have real parts identically zero and the 
imaginary parts tend to the same value. On the other hand, the same solution 
is obtained [llj by taking the limit of the breather solution when a goes to zero 
(see section 2 below for the explicit expressions) . In this case the poles have an 
identical imaginary part while the real part is changing. 

The main purpose of this paper is to see that these numerical instabilities 
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do occur and even regular numerical methods as pseudospectral and finite dif- 
ferences fail when the time simulations become long enough. In the case of 
the finite differences, studied in section 4, two spatial discretizations have been 
used, with different discrete invariants associated to each. The motivation is 
to highlight that the numerical approximation to the two pole initial condition 
separates from the right solution and chooses one of the two possible branches, 
either two independent solitons or a breather. It is interesting to observe that 
the pseudospectral method, presented in section 3, is the most efficient studied 
here. It captures the logarithmic separation of the two humps for much longer 
times than the finite difference schemes. Nevertheless if a choice of a low num- 
ber of harmonics or a too big step in time is made, then a wrong behavior also 
appears in the case of pseudospectral methods. 

In the literature it is found that the accuracy of the numerical solutions 
is directly related with the consistency order of the numerical approach. For 
example the dynamical instability of the breather solutions of this equation as 
well as the stability of the double pole was numerically studied for reasonable 
long time intervals in Ref. [17] . The convergence of the methods only guarantee 
good approximation in short time simulations. But these limitations may be 
overcome by fixing some invariant quantities, which are constant for the exact 
solution, along the evolution in time. These quantities are the discrete equivalent 
to the conservation laws of the continuous equation as the mean, the norm 
or the energy. 

There are several works concerning the stability of a soliton with one hump. 
Also rather progress has been done in the case of multisolitons when the humps 
are widely separated from each other, and even more recently about the collision 
of two solitons, one fast and narrow and the other one slow and broad. However 
very little is known in the case we are considering in this paper. The available 
analytical techniques do not seem to apply in this case |18j . We want to empha- 
size that the instability observed in this work do not lie in the arbitrary growth 
of the error, but in the incorrect behavior of the numerical solutions that fall 
within the orbits corresponding to different solutions in the phase space. In this 
paper we give numerical evidence that supports the difficulty of the problem. 

2 Initial condition 

We will study numerically some particular vanishing at infinity solutions of the 
modified Korteweg-de Vries equation (JlJ in 1+1 dimension, (x,t) 6 [0, T] x M, 

lim u(x,t) = 0. (4) 

We are interested in a two parameter family of exact solutions that were 
first obtained by Wadati in Ref. [11]. They can be written as 



u(x, t) = 2p- sech (j8 (x + -ft)) x 
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cos[$(x,i)] - I - ) sin ($(x,i)) tanh (^(1 + 7*)) 



1 + ( sin 2 t)) sech 2 (j8(a; + 7*)) 



(5) 



with 7 = 3a 2 - /3 2 , <S = a 2 - 3/3 2 and $(a:,t) = a(x + St) - taa _1 (/3/a). They 
can also be written as 

^,t) = 2^tan-^ ffil^l§ ) (6) 
\F(x,t,a,p)J 

with the auxiliary functions -F(a;, t, a, /3) and G(x, t, a, /?) given by 

F(x, t, a, 0) = cosh (/3 (3to 2 - t/3 2 + x) ) , 

/3sinfa(x + t(a 2 -3/3 2 )) -tan" 1 f^)) (7) 

G(x,t,a,0) = ^- ^L. 

a 

Their shapes correspond to a breather with a group velocity —7 — /3 2 — 3a 2 
when a^O (Figure QJi) that degenerates to a double pole solution when a — 
(Figure [TJj) . This second case is interesting because it behaves asymptotically 
as a pair of independent solitons of the mKdV Eq. (Q}, one facing up and the 
other one facing down, 

. . 4^ (cosh (/3a; - /3 3 t) + /3 (3fi 2 t - x) sinh (/3x - /3 3 t)) 
ii(x t) — —— (8) 

2 (/3x - 3/3H) 2 + cosh (2(3x - 2/3H) + 1 

This double pole solution has three local extremes with one of them decaying 
to asymptotically in time. The other two extremes separate with a distance 
l{t) that becomes logarithmically large in time when t is big enough, see Ref. 

EL 

„ . 21ogf4^ 3 -t] 

«(*)»— l , (9) 

This result is easily found by taking into account the point-wise convergence of 
the function ([SJ to when time goes to 00, and checking the position where the 
derivative of the exponential dominant terms vanishes. The property ^ will 
be very useful to check the accuracy of the numerical methods for t big. 

One of the most important properties of this type of partial differential 
equations is the existence of infinite conservation laws. For example the solution 
u of the mKdV equation on the real axis x G = R preserves the mean in space 
I(u), the mass ||u||2 and the energy E(u) over time, as detailed below. These 
statements are straightforward by integration by parts as well as taking into 
account the boundary conditions ([3]), 



I(u) = / u(x, t) dx, 
Jn 

\\u\\ 2 = / u 2 (x,i)dx, (10) 
Jn 

E(u) = / [u 4 {x,t) -u 2 x {x,t)] dx. 
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In this paper we will show by computational experiments that in long time sim- 
ulations the numerical schemes that preserve exactly some of these conservation 
laws, Eq. dTU|) . are much more robust than the schemes that only take into 
account the consistency order of the method. The numerical data have been 
compared with the exact analytical expression for the travelling wave solutions, 
written by ([5]) and ©, in order to measure its accuracy. In practice the ini- 
tial condition provided by the experimental measurements, which have to be 
investigated numerically, may differ from those functions. The use of a powerful 
numerical scheme plays a crucial rule in long time experiments. 

Due to the fact that the double pole initial condition © is very regular and 
no oscillations are present, both type of methods seem to be appropriate to 
simulate its evolution. On the other hand, when the highly oscillating initial 
condition Eq. (|5|) is considered, we postulate that the pseudospectral method 
will be the most powerful because an acceptable number of harmonics can rep- 
resent with high accuracy the shape of the solution. In this case the finite 
difference methods need a too high number of nodes to generate an acceptable 
approximation of the solution, and the global error of the numerical approxi- 
mation will increase exponentially due to the big Lipschitz constant induced by 
the oscillations of the pulse. 



3 pseudospectral method 

As a general rule, the pseudospectral collocation methods, Refs. [15], [3DJ and 
[21] . are very suitable to approximate the travelling wave and breather type 
solutions of a partial differential equation. The oscillating profile of a function 
and its derivatives can be reproduced quite accurately by a manageable set of 
harmonics Ref. [32]. Due to the limitations of the practical discretization in 
space, the real axis R has been substituted by a long finite interval f2 = [-L, L] 
(the numerical simulations have been made for L = 40) and the boundary 
conditions ([4]) by periodic ones, 

u(-L) = u{L), u x (-L) = u x (L). (11) 

These conditions suggest to use an orthogonal basis made up of the complex 
exponential periodic functions, {4>j(x) = e lUjX }, with Wj — 2irj/(2L), j 6 Z. 
The Fourier series expansion of a function u(x,t) £ L 2 is defined by 

oo \ f L 

u(x,t) = Uj(t)(f>j(x), where Uj(t) = — / u(x, t)<f)j{x) dx. (12) 

j=-oo J - L 

At initial time, t — 0, the infinite expression (|12|) is replaced by some trun- 
cated series at the collocation points, x n = 2nL/N, n = —N/2, . . . , N/2 — 1, 
to adapt it to the discrete numerical scheme, 

N/2-1 N/2-1 
U^=U (x n )= ]T U^Mx n ), Uf = ^ £ MXn)h(x n ) (13) 

j=-N/2 n=-N/2 
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Fig. 2: Theoretical separation l{t) of the humps of the double pole solution for 
(3 = 1 (solid line) compared with the result of pseudospectral approx- 
imation for N = 2 8 (dots), N = 2 9 (triangles) and N = 2 10 (circles) 
number of collocation points and time step At — 10 -3 . 

The advantage of this approximation is that it can be notoriously acceler- 
ated by the well known FFT algorithm when N — 2 9 , gGN, collocation points 
are taken. In our work we have ran simulations for N = 2 8 , N = 2 9 and 
N = 2 10 points to check the precision of the discretization and the robustness 
of the method. These calculations have been computed by a FORTRAN sub- 
routine provided by Netlib repository [23] . Some authors (24] [25l [26] separates 
the linear and the nonlinear parts of the right hand side of the equation (p} and 
later on they introduce an exponential integrant factor in the solution to deal 
with the linear part. The advantage of this reduction is that the only term to 
be discretized by Fourier transform or by finite differences is the nonlinear one 
and its norm has order O(N) instead of 0(N 3 ). In this context the mentioned 
strategy would not produce significant improvement in the results because the 
wave-number of the functions of the orthogonal basis used to approximate the 
solutions is not necessary high. The situation is radically different when solu- 
tions that involve high frequency radiations are investigated, then it is really 
interesting to avoid the computation of the derivatives of the linear part by 
using the mentioned integrant factor technique. 

Let us define = (u^^ 2 , £^vy 2 -i) as tne vector that stores the 

initial data, and « u{tk) as the successive approximations of the solution 
at time ifc. Then, the discrete Fourier transform can be considered as a linear 
operator t/( fc ) = represented by a Vandermonde N x N matrix T with 

components Tj, n = 0( i "^" 1 X n "^ ^) and 9 = e -^/ N . 

Taking into account the orthogonality (<fij, <pi) = Sij of the periodic functions 
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Fig. 3: The same as Fig [2] for f3 — 1, N — 2 9 and different time step sizes, 
At = 10 -2 (dots), At = 1(T 3 (triangles) and At = 1(T 4 (circles). 



4>j in L2[— 7r, 7r], and the expressions of their space derivatives, 4>'j{x) — iujj4>j(x), 
we substitute the expansion (|12p in the solutions of Eq. ((TJ) obtaining an ODE 
system for the Fourier coefficients, 

d t u 3 {t) = i [ujujit) - 2w j t) J -(t)] , j = -f , f - 1, n , 

u(0) = Fu(x,Q). v ; 

Here v = J-v is the Fourier transform of the nonlinear term v — u 3 . Substituting 
in Eq. (|14j) the continuous solution u and its power v — u 3 respectively by the 
iV-dimensional vectors U and V with the values at the collocation points, we 
obtain the equivalent matrix expression for (|14l) . 

d t TV = - [V 3 TU + 2VTV] . (15) 

The matrix T> is diagonal with components djj — iujj and represents the spatial 
differentiation. Considering the linear antisymmetric operator J = T~ lr DT ', 
where J — — J , the previous relation (|15l) corresponds to the motion equation 
of the Hamiltonian defined by 



(16) 



Taking advantage of the conservation of the linear and the quadratic invariants 
Eq. (|10l) by the symplectic integrators, as far as the consistency order of the 
scheme, we suggest to use the implicit midpoint rule for discretization in time, 
together with the pseudospectral method for discretization in space. We postu- 
late that this property about the discrete method will lead to an improvement of 




3 pseudospectral method 







v(t) 

74.5 r 



74 



73.5 



73 



CQOOOOOOOOOOOOOOO°OOOOOOOOOOOOOOOoOOOOOOOOOOOOC 

""iililiiiilllli 



7V=2 8 
7V=2 9 



7V=2 



10 



100 



200 



300 



400 



500 



Fig. 4: The group velocity of the approximated breather (a = 5 and /3 = 1) 
for At = 1CT 3 and different number of points N = 2 s (dots), TV = 2 9 
(triangles) and TV = 2 10 (circles), versus v = 3a 2 — (3 2 = 74. 



the time interval where the numerical solution approaches the exact one. Even 
though the midpoint rule is not a method with a high order of consistency, the 
conclusions obtained from its results will be representative and clear for the 
purpose of this paper. In addition the convergence for each time step is easily 
fulfilled, contrary to the family of Newton's type methods that are much faster 
than the midpoint rule but they often require deep analysis in order to avoid 
the problems caused by local minima. 

Let be the discrete Fourier transform of U^ k \ then the midpoint rule 
applied to the pseudospectral discretization of Eq. (fTSj) reads, 

jj( k+ i) = jj(k) _ At r v s(j(k+i) + 2V y( k+ ^ ; (17) 

where u( k+ ^) = 0W + f7( fe+1 ))/2 and ^O+i) = (y(*0 + t>( fe + 1 ))/2 are the 
interpolated vectors. Putting together the terms t/( fe+1 ) at the left hand side 
of Eq. pT|) the norm of the subsequent functional associated to the fix-point 
scheme is considerably reduced. A successful implementation is provided for 
time steps of size At « 10 -3 , where an acceptably cheap computational cost is 
involved. We get the following iteration formula 



I - — V 3 ) U {k) - 2Afflv( k+ ?) 
2 



(18) 



Notice that the matrix [I — (Ai/2)Z) 3 ] is diagonal and its inversion is directly 
made. The real and the imaginary part of the system (|18|) are separated in 
different equations for the practical implementation. 
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Fig. 5: The group velocity of the approximated breather (a — 5 and j3 = 1) 
for N = 2 9 and different time step sizes At = 10~ 2 (dots), At = 10~ 3 
(triangles) and At = 10~ 4 (circles), versus v = 3a 2 — (3 2 = 74 . 



In the numerical simulations with the double pole initial condition we have 
found a very long regular behavior of the evolution. The Figures [2] and [3] show 
the logarithmical evolution in time of the separation distance between the main 
humps of the approximate solutions ((9|) with respect to the exact l(t). The 
increment of the number of points from N = 2 9 to N — 2 10 does not improve 
the accuracy as much as the reduction of the step in time from At = 10~ 3 to 
At = 10~ 4 which provides a good agreement with the exact solution. Only 
when the number of harmonics is too low, N ~ 2 8 , then the double pole initial 
condition can break on two independent solitons Eq. ([2]) as it is observed in 
the doted trajectory of Figure O which is almost linear and corresponds to a 
successive separation of the two main humps of the solution with different but 
nearly constant velocities. 

The results for breather initial conditions with a = 5 as internal oscillation 
frequency have been summarized in Figures [4] and [5] The choice of relaxed 
restrictions on the time and space steps At > 10~ 3 or N < 2 9 breaks up on 
the damage of the numerical solutions and the corresponding approximation 
to the conservation laws of the continuum equation (jTUJ) . However a not so 
expensive conditions for the scheme (|18p as At w 10~ 4 or TV = 2 9 guarantee 
a very accurate behavior and preservation of (|10[) during long time intervals, 
t € [0, 500], as it can be observed in Figs. 0]and[5] 

It is important to explain carefully the results captured on the figures ((4]) and 
f5"| in order to avoid confusions on the interpretation. In the first one different 
number of collocation points have been considered while in the second one it is 
investigated the accuracy of the solutions for different time step sizes. In both 
cases the decreasing of the number of points of the mesh leads to the same effect, 



4 Finite difference methods and some discrete invariants 



11 



which is the faster loss of the convergence to the exact solution because of the 
jump to another family of solutions. 

4 Finite difference methods and some discrete invariants 

A wide number of papers have appeared in the last years investigating the 
accuracy of finite difference methods Refs. [37], [8], [S] applied to nonlinear 
dispersive partial differential equations and comparing them with other type of 
discrete schemes as pseudospectral or Adomian decomposition. The goal of our 
analysis is to give evidences of the connection between the conservation of some 
invariants of the discrete scheme and the improvement on the convergence of 
the numerical method to the exact solution, in order to avoid the instabilities 
we mentioned in the introduction. 

In order to implement a numerical method which approximates the solu- 
tion of the initial value problem ([1]) we first define a spatial discrete grid of 
points x n = —L + nAx, n = 1,...,N, and the corresponding vector with 
the approximation to the solution in that points, U(t) — (J7x(t), . . . , Urf(t)), 
where U n (t) ~ u(t,x n ). Let's be M(u) = — Su — B(u) the differential op- 
erator that concentrates the spatial derivatives of the mKdV equation ([1]), 
where the linear part is assumed by S(u) = d x (u) an d the nonlinear part by 
B(tt) = 3ud x (u 2 ) = 2d x (u 3 ), there are several choices for an operator in finite 
differences Mn(-) = — Sn(-) — Bn(-) that approximates M(-) with different 
order of consistency. The time evolution of the components of the vector U is 
governed by the following system of ODE, 

d t U = M n (t,U), (19) 

provided by the periodic boundary conditions (jlll) in [— L, L], that can be writ- 
ten as 

U n+jN =U n , Vj G Z, < n < N. (20) 

The accuracy of the numerical approximation U with respect to u depends 
on the consistency order of the operator Mn(-)i which consists of the sum of 
some suitable band matrices. Let's define first the forward difference matrix D^, 
where only nonzero components in each row are di_i+i = 1/Ax = —di t i, i = 
1,...,N, and djv,i = 1/Ax in the corner due to periodic boundary conditions. 
It represents an approximation of the spatial derivative, d x u w D^[7, of the 
first order O(Ax). Using the expression of the matrix we can define the 
classical discrete approximations of the successive derivatives: 

a) Backward differences, = — (D^) T , which approximates d x with order 
0(Ax). 

b) Central differences, D£ = (D^ + )/2, which approximates d x with order 
0(A 2 x). 

c) Central differences, D£ = • , which approximates d x with order 
0{A 2 x). 
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d) Central differences, D| = DJ-D^, which approximates 3 3 with order 0(A 2 x). 

It is straightforward to check that and D£ are antisymmetric matrices 
and D2 is a symmetric matrix. This observation will be important to guarantee 
the conservation of some invariants of the discrete schemes along the time. Now 
we write the linear part as 

S N £7 = D|£7. (21) 
However the nonlinear part Bn((7) admits some suitable alternatives as 

B^{U) = 2D5t/ 3 or Bjf(U) = 3UiD5t/ 2 . (22) 

Here the vectors U 2 and U 3 are defined respectively by U 2 = (Uf,..., Uff) and 
U 3 = (U 3 , . . . , Upf) and the matrix Ui is a diagonal matrix with the components 
uu = Ui. From now on we will suppress the t dependence in the notation of 
Mn due to the autonomous structure of the mKdV equation. 

Next we will choose a numerical scheme for the discretization in time, where 
Un ~ u(x n ,tk) approximates the exact solution at time As we have 
mentioned in the previous section, the use of a symplectic integrator in the 
forward time step is suitable for this kind of evolution equations. They pre- 
serve the linear and quadratic invariants of the continuous equation, as far as 
the consistency order of the numerical method (Ref. [28], [[29]] ). Because of 
that we use the implicit midpoint rule that evaluates the flow of the equation 
M5j([/) = -S N [/-B^([/), r = 1,2, at the vector lf( k+ ^ = (U^ + U^ k+1 '>)/2 
placed between the approximations at tk and tk+i, 

v {k+i) =U (k)_ At f SNU (k+i) + B^(u( k+ ^, r = l,2, (23) 

These are finally the two implicit numerical schemes in finite differences to be 
analyzed. 

4.1 Convergency analysis 

The finite difference schemes proposed above will be C-stable (Ref. 28 ) if the 
logarithmic norm of the Jacobians, //2 (M£/([/)) , r = 1, 2, are bounded inde- 
pendently of the grid spacing of the evolution operators for every U belonging 
to the segment between C/ fc+ 2 and u(x, t + At/2). In both cases considered in 

Eq. (|231) this norm is bounded by [i 2 (M|/(Z7)) < max n (Un + ^) 2 . In Ref. [S] 
it was proven that in the continuous case this is true for initial conditions U° 
in Loo. 

To prove the C-stability of this implicit method we consider two different 
set of numerical approximations {U^}k>o and {U^}k>o, associated to two 
slightly different initial conditions, and a bound for the squares of the numerical 

solutions \{U,{ k+ ^) 2 \ < n, for all k > 0. Then 
^(fc+i) _ fj(k+i) = uW _ (j(k) + At ^ (u(k+i)^ _ M r^ f(j(k+m\ . (24) 
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Now multiplying the expression (1241) by and applying the 

mean value theorem to one gets that 

- ||2 = i||C/W - J7«|| 2 + 

At {M^ (u( k +i)} - M r N (u( k +^ (u( k+ i) - ll( k+ ^ < (25) 
i\\UW - UW\\ 2 + Atfi\\u( k+ i) - u( k+ i)f. 

Using the triangular and the Cauchy-Schwarz inequalities we deduce 

|| C/ (fe+i)_ ? y(/c+i)|| < ^jj(k) ^-(fe)||2 ! Af/z (njj(k+i) _ [/(fc+i)]! + _ {>( fc )||^ 2 

(26) 

To end the proof for the C-stability we consider the relation 

H^fc+i) _ £(fc+i)|| ^(fc) _ ^(fe)|| < m^n^fc+i) _ f/^+Df, _ ^«|| 2 }, 

(27) 

and now we can take common factor of some terms to conclude that 



H^fc+i) _ ,| < sj l+^ 2 \\UW - U^l < 3 M A</2 < 1. (28) 

The consistency of the full scheme after one step in time is proven following 
Ref. [28] defining (j( k + 1 ) = (ji k + 1 ) + AtM^U^+i)) as the numerical solution 
at time t k +i being U^> — u{tk) the exact solution. Then the global error is 

(3(t k+1 ) = - u(tfc+i) (29) 

By defining the truncation errors corresponding to the midpoint rule and the 
trapezoidal rule respectively as 



d x (t k+1 ) = u(t k ) - u(t k+1 ) + AtM 



N 



u{t k ) + u{t k+1 ) 



d 2 (t k +i) = u(t k ) - u(t k+1 ) + ^ (u(t k )) + M T N (u(t k+1 )) 



(30) 



then the global truncation error can be written in the following manner 
P(t k+1 ) = dx(t k+1 ) + At 



M r(^)^ {k+1) )_ M rJ^ k ) + u(t k+1 ) 



2 

'(31) 

Now using again the mean value theorem with the operator MJ^, multiplying 
Eq. (JSIl) by J/( fe+1 ) - u(t k+1 ) and using the Holder's inequality we obtain 

\\U^ - u(t k+1 )\\ 2 < (l-^) l|di(**+i)lla, A*At<2. (32) 
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To complete this argument it can be deduced by Taylor expansion that, for 
sufficiently smooth u, the function d\ can be bounded as ||<ii(£fc+i)||2 = 0(A 3 t+ 
AtA 2 x) and consequently for < fiAt < 1 

ll/3(t ^+ l)l12 = 0((At) 2 + (Ax) 2 ). (33) 

This result finally gives us the proof of the convergence of these schemes for a 
given time interval [0, T] whenever the steps Aa; and At are chosen small enough 
depending on T. But, as we will see later, this statement does not prevent the 
method to corrupt the solutions for longer time intervals. 



4.2 Discrete invariants 

In the literature there are many works that emphasize the importance of design- 
ing numerical methods that preserve the invariants related to the conservation 
laws of the continuous model. For this purpose different strategics as geomet- 
ric integration [3UJ or projection methods [3T] have been developed. In this 
section we investigate the transcendence of preserving some invariants by the 
structure of the full discretization schemes (|23[) for improving the stability of 
the numerical solution corresponding to the initial conditions, either (J5]) or Q. 
We will focus our analysis on the conservation of the following quantities, 

£ 1 {U) = AxJ2u n ~I(u) 1 

n 

d(U) = Ax^(U n ) 2 = (U,U) k \\uf, 

n 

C 3 (U) = AxY, [(Un) 4 - (U n - £4-i) 2 ] = {U 2 , U 2 ) - <D+C/,D+[/) « E(u). 

n 

(34) 

These expressions are the discrete approximations to the continuum invariants 
(fT0|) of the mKdV equation, by the discretizations (|23|) presented in the pre- 
vious section. The previously proposed discrete schemes present interesting 
advantages with respect to the conservation of the discrete invariants (|34[) . The 
first scheme, U {k+l) = U (k) + AtM^(u( k+ ^), preserves the linear quantity 
d{U k ), Mk > 0, 

& = Ci (t/ (fe) ) + At d (mJj (V( fe+ 2))) = C^U^). (35) 

Here we have used that the sum of the components of the columns of the 
matrixes Sn and Bjj involved in the definition of is 0, consequently 

£i(Mjj([/( fc+ 2 ))) = 0. This property is not satisfied by the second scheme 
jj(k+i) = jj(fe) + AtM^(t/( fc+ ^), however this discretization preserves C 2 {U k ), 
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Vfc > 0, as it can be easily shown, 

c 2 {u (k+1) ) = {u ik+1) ,u (k+1) ) = {u (k) + u (k+1) - u (k \u( k) + u {k+1) - u {k) ) 

= £ 2 (C/ (fe) ) + '- U^, + C/( fe+1 )) + <!/(*), C7-(*+i5 - [/(*)) 

2 \ (36) 

The last simplification arises from the skew-symmetry of the matrixes Sn and 

B|j({7), which are involved in the definition of M^j. 

The linear and quadratic invariants shown above is all we can expect to 

be preserved by a symplectic method along the time steps. Their numerical 

precision will be of the same order as the tolerance imposed on the numerical 

solution of the nonlinear system (|23l) . In addition we can prove that the exact 

time integration of the spatial discretization dtU = Mjj({7) would maintain 

constant along the time the invariant C3QJ), 

d t C 3 (U) = d t ((U 2 , U 2 ) - (D+U, B+U)) = 4 (d t U, U 3 ) - 2 (B+d t U, B+U) 
= 4 (d t U, U 3 ) + 2 <(ftC/ T Di ), B+U) = 2 (d t U, (2U 3 + B C 2 U)) 
= -2 (BlU + 2BIU 3 , (2U 3 + B$U)) 

= -2 (2(D§[/,C/ 3 ) + (D%U,B%U) +4(DJ U 3 ,U 3 ) +2 (D^f/ 3 , D°C/)) 
= -2 h (B%U, U 3 ) - (u, BZ T B$U^ + 4 (D^ 3 , U 3 ) + -2 {U 3 , B\B%U)\ = 

(37) 

The fact that B\ and D^D^ are skew-symmetric matrixes and B\ T B 2 — 
— D^Dj = — D3 has been used to proceed with the cancelations in (|37|) . The 
approximation to the first spatial derivative in £3 (£7) by B + and the definition 
of Mjj([/) are crucial for this purpose. However it is easy to prove that the 
scheme dtU — M^j(J7) will not reach the condition (|37|) . 

Summarizing, the spatial discretization d t U = Mjj([7) preserves the invari- 
ants Ci(U) and C^(U). The first one of them is also preserved by the full dis- 
cretization in time by the midpoint rule. On the other hand the full space-time 
discretization corresponding to M^([7) preserves £2(1/). The consequences of 
these results on the stability of the numerical solutions of the type (0) on long 
time intervals will be discussed in the next section. 



4.3 Implementation of the numerical scheme 

In this section we will describe the implementation of the schemes (|23|) that 
give approximations to the double pole solutions ([8]) of the mKdV equation (fTJ). 
The results will be compared with the exact solution in order to measure the 
efficiency of the numerical methods and to obtain conclusions. 

The implicitness of the schemes makes necessary to apply them together 
with a numerical method for solving systems of nonlinear equations as a fix- 
point iteration Ref. [32 . First we rewrite the general method (|23|) as a sum of 
its linear and nonlinear part, 

C/ (fc+1) = U^+AtMl, ([/( fe+ ^)) - U™-At [s N C/( fe +^) + B r (V^))] , r 

(38) 
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Fig. 6: (a) Finite difference approximation by Mj^ of the double pole (/3 = 1) 
for At = 4 ■ 10~ 2 and Ax — 10 _1 (dots) versus an exact breather of the 
mKdV for a — 0.0982 and (3 — 1 (continuous line), (b) The same for 
At = 1 ■ 10~ 2 versus two exact solitons of the mKdV with amplitudes 
v = 0.95 and v = 1.05. 



The time step At in (|38|) should be chosen small enough to assure that the 
operator AiMJ sr ( / J fc+ 2 ) is contractive with respect to U^ k+1 \ In practice the 
grid size Ax = 1/N imposes a very restrictive choice of At to guarantee that 
the norm of ||A£Sn|| = 0(N 3 At) is lower than 1. This circumstance suggests 
to treat the implicit term u xxx replacing scheme (|38|) by the following formula, 
as was made in Ref 1331. 



/j(fe+i) 



I 



At. 



I 



At. 



s N i u {k) - Am 



fe+i) 



(39) 



Now the norm of the fix-point operator is independent on N. 

When a double pole initial condition ((SJ is chosen, the results of the simula- 
tions fall into a very reach casuistic, and gives a strong hint about the instable 
nature of this type of solutions in the continuum case. As it was mentioned in 
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the introduction, this particular solution is related to the fact that the reflec- 
tion coefficient that appears in the inverse scattering formulation has a second 
order pole in the imaginary axis, Ref. [TT]. This situation can be understood 
as a limiting case of solutions either made as superposition of two independent 
solitons (the poles are different and are locate on the imaginary axis, j3\i sa {3 2 i) 
or as a breather of low frequency (the poles are different and symmetric with 
respect to the imaginary axis, ±a + /3i, a << 1). 




Fig. 7: Evolution of the solution by finite differences (Ax = 10 _1 ) for differ- 
ent At and discretization (left) or (right), from a double pole 
(dark regions) to two independent solitons (gray regions) or a low oscilla- 
tion breather (dashed region). The big dots indicates when £3 abruptly 
oscillates. 

This complex scenario, where structurally different type of solutions are 
distinguished by the slight change on the choice of the parameters a and /3, is 
captured by the behavior of the approximations produced by finite difference 
methods. We observed the following phenomena: 

a) For long time experiments, the double pole solution is better approximated 
by the scheme where the spatial discretization is approached by the operator 
M^j than the other choice (dark region in Figure [7]). A time step as small as 
At p» 0.02 has to be taken to guarantee stability. 

b) After reaching a transient time, the initial double pole shape is trapped 
by the orbits of two independent solitons when spatial discretization is 
used or Mf^ with too small time step (gray region in Figure [7]). An example 
of this phase change is shown in Figure [6](b) , where Mf^ spatial discretization 
with At = 10~ 2 , causes that the initial double pole shape is captured by two 
independent solitons with amplitudes (3 = 0.95 and (3 — 1.05. A similar behavior 
would have been obtained for the discretization Mjj with any time step after a 
prudential time interval T > 50. 

c) When M^j spatial discretization is used together with not so small time 
step (At > 0.02), then the original double pole solution, after going through the 
shape of two independent solitons for a short period of time, leads to a breather 
solution with a small a parameter (dashed region in Figure [7]). This jump on 
the shape of the solution is detected at the same time by a breakdown of the 
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energy Cs(U) as is shown by the position of the big dots in the right diagram of 
Figure [71 The use of discretization avoids this pathological behavior due 
to the intrinsic preservation of the energy. An example of the trapping by the 
orbit of the solution ([5]) with parameter a = 0.0982 is shown in Figure UJa). 

The conclusion of this section suggests that the methods that conserve the 
mean L\ and the energy £ 3 preserve the numerical double pole solution from 
jumping to the orbit of a breather in the phase space of solutions. This patho- 
logical behavior is specially caused by the use of M^j discretization together 
with a roughly choice of time step At. On the other hand the discretization 
for M^j guarantees good results for a longer interval of time than when a 
small enough At is chosen. It means that the conservation of the discrete £2 
is an important feature to be taken into account by a numerical method that 
reproduces the exact behavior of the solution for long time simulations. 

5 Conclusion 

In this paper we have studied the response of the pseudospectral methods and 
the finite difference methods when they are applied to a particular family of 
solutions of the mKdV equation. These solutions, known as breathers, depend 
on two parameters, a and /3, which control respectively the internal oscillations 
as well as the effective group velocity of the waves. When the a parameter 
vanishes, then a special type of solution called "double pole" appears. Its shape 
is characterized by two main humps, one up and the other one down, which 
progressively separates with a logarithmic velocity. This solution can be con- 
sidered as a limiting case of the breather solution when a — > and also it can 
be seen as the superposition of a soliton and an antisoliton. Therefore from the 
the point of view of the inverse scattering theory the double pole solution has to 
be seen as a degenerate case where three type of solutions are very close. This 
circumstance explains that when perturbations are considered the correspond- 
ing solutions are captured by any of the three possible orbits. The first two 
branches are either to follow a breather with a > or a soliton type solution 
(a = 0). In this second possibility another two orbits are possible, either the 
wrong branch where the amplitudes of the two humps become different and they 
separate by a linear speed, or the right one where the amplitudes tend to the 
same value and the separation is at a logarithmic rate. 

The type of smooth solutions of the mKdV equation considered here clearly 
suggests that pseudospectral methods are the most appropriate for describing 
their dynamics. In fact, the functions of the basis of the Fourier discrete ex- 
pansion adapt themselves very easily to the oscillations of the solution, and it 
is not necessary a high amount of harmonic functions to obtain an accurate 
approximation. However, this intrinsic oscillations of the breathers involve a 
big Lipschitz constant even for a > not necessarily large, and they make the 
finite difference schemes not feasible to obtain good approximations in this case. 
Therefore we have only used them for the simulations that have the double pole 
as initial condition. 



6 Acknowledgment 



19 



It is well known that the mKdV equation has an infinite number of conserva- 
tion laws. As a consequence of it, we propose numerical schemes that preserve 
some of them with the aim of getting accurate approximations to the solu- 
tions in long time intervals. In order to increase the possibilities of success, the 
time advancing is reached by a symplectic method as the midpoint rule. This 
implicit method involves a system of nonlinear equations that is conveniently 
manipulated to be economically solved by using a fix-point procedure. 

The results of the simulations made by finite difference methods bring to 
light the relation between each family of solutions and the preservation of some 
of the conservation laws. Keeping constant the discrete invariants equivalent to 
Ci{U) and specially to the energy C^{U). The discretization that keeps constant 
the equivalent to the integral of the function prevents the initial condition to 
degenerate from a double pole to a breather. Unfortunately after a short time 
the solution drops in two independent soliton-antisoliton shape. However the 
discretization that keeps constant the equivalent to the integral of the square 
of the function, £2(1/), reproduces for a rather longer time the double pole 
solution with a good accuracy. Nevertheless the approximation deteriorates for 
later times and jumps either to a breather type solution when the time step is 
rough At > 2 • 10~ 2 , or to a soliton-antisoliton solution when the step is small 
enough. 

The pseudospectral method has been successfully developed with both types 
of initial conditions, the double pole as well as the breather. Taking At = 10~ 4 
and N = 2 9 for a double pole, and At = 10~ 5 and N = 2 9 for a breather, 
the method reproduces with extremely good efficiency the exact solution dur- 
ing long intervals of time, t e [0,500]. An important remark is that in these 
conditions, the pseudospectral discretization keeps constant the three discrete 
quantities investigated in this work, as far as the precision order of the method 
allows. Considering implementations of comparable computational time for the 
pseudospectral and finite difference methods, the first ones return good results 
for time intervals of one order of magnitude longer than the second ones. 
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